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Abstract. In this paper we propose a method to couple two or more explicit numerical schemes 
approximating the same time-dependent PDE, aiming at creating new schemes which inherit advantages 
of the original ones. We consider both advection equations and nonlinear conservation laws. By 
coupling a macroscopic (Eulerian) scheme with a microscopic (Lagrangian) scheme, we get a new kind 
of multiscale numerical method. 
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1. Introduction 

In this paper we propose and analyse a method to couple two or more numerical schemes which approximate 
the same time-dependent PDE. The result is a new scheme which inherits advantages (and drawbacks) of the 
original ones. In what follows, we will discuss how the combination might allow to improve each of the schemes, 
and try to derive practical coupling strategies. Roughly speaking, we combine the existing schemes by a convex 
combination in such a way that each scheme is influenced by the others. In the end, we get as many approximate 
solutions as the number of the considered schemes, but all of them are modified by the combination. In the 
basic version of the algorithm the coupling is performed everywhere and at any time. Nevertheless, local (in 
space or in time) couplings are also possible. 

For illustrative purposes, we focus on the coupling of two schemes only. We consider both the case of two 
macroscopic (Eulerian) schemes and, more interesting, the case of a microscopic (Lagrangian, particle-based) 
scheme blended with a macroscopic scheme. In the latter case we get a new kind of multiscale numerical method. 
As a test problem, we consider the following advection equation in conservative form 

( ut(x,t) + (A(x) u(x,t)) x =0, x £ t > 0, 

( u(x, 0) = u(x), x £ 
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where u : R d —> R is the initial datum and A : R d —> R d is a bounded and Lipschitz continuous vector field. 
This equation describes the transport of a certain amount of mass in R d , drifted with velocity A. The function 
u represents the density of such a mass, so that the mass in any set E C R rf at any time t is given by 

M{E,t) = / u(x,t ) dx. 

J E 

In order to introduce a multiscale framework, it is useful to recast the problem at microscopic scale. Denoting 
by V £ R d the position of a generic mass particle, the Lagrangian counterpart of ([I]) is the ODE 


V = A(V). 


( 2 ) 


with suitable initial condition 'P(O). 

The extension to the one-dimensional nonlinear conservation law 


u t (x, t) + f(u(x, t)) x = 0, x £ R, t > 0, 
u(x,0)=u(x), x £ R, 


( 3 ) 


will be also discussed. 

The coupling idea investigated in this paper can be pursued whenever one has at his disposal two or more 
numerical schemes with complementary advantages and drawbacks , resulting in a blended scheme which shows a 
combination of these features. The idea is mainly inspired by the multiscale technique introduced in j4], which 
can be applied whenever one has two equations modeling the same physical phenomenon at different scales, 
and the velocity field explicitly appears in both equations. The blend is obtained by a convex combination of 
the two velocity fields, an then both the macroscopic density and the microscopic particles are transported by 
means of the new common multiscale velocity field, which, in turn, depends on both density and positions of 
the particles. Along the same line, the paper 3 blends the Brownian motion with the heat equation. Note that 
in these papers the goal is to improve the description of some physical phenomenon, by catching at large scale 
some effects ultimately triggered by small scale features, otherwise invisible. Conversely, in the present paper 
we blend directly the solutions of the schemes rather than the velocity fields, with the goal of improving the 
quality of the numerical approximation. 


Relevant literature 

The literature about numerical schemes for the advection equation, including multiscale numerical methods, 
is huge and a complete review of it is out of the scope of the paper. We will refer the interested reader to 
the book 13 for an introduction to the equations ([!]) and ([3]), and the basic Eulerian numerical schemes to 
solve them. We also mention the book [5] for a survey on the best known techniques to couple microscopic and 
macroscopic descriptions of physical phenomena. In the following we only discuss the schemes that can be put 
in relation with the one proposed in this paper. 

• The proposed scheme follows the same philosophy of Crank-Nicolson method and 0-methods in general. 
However, there are two important differences: first, our method is explicit and combines two explicit 
schemes, while 0-methods always combine the explicit and the implicit version of the same scheme. 
Second, our method has two free parameters instead of one. 


Filtered (hybrid) schemes [8} 9 17 suitably combine a monotone scheme with a high-order scheme, in 
order to get a high-order approximation minimizing spurious oscillations. In this case the coupling is 
restricted to Eulerian schemes. 

The Particle-in-Cell scheme jl , often used in plasma physics, combines a Lagrangian scheme with an 
Eulerian grid. Roughly speaking, Lagrangian particles are independently tracked, but an underlying 
grid is used to recover averaged quantities. At each time step, those quantities are used for additional 
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computations which finally give the force acting on each particle. Then the procedure is repeated for 
the next time step. 

The Particle Level-Set (PLS) method [6] is a modification of the standard level-set method 18 for 
interfaces tracking, introduced to efficiently reduce the numerical diffusion appearing at regions of the 
evolving front with high curvature. Lagrangian markers are posed in a neighborhood of the interested 
region, and moved independently by means of a high-order ODE solver, according to the characteristic 
flow of the level-set equation. Whenever a particle escapes from the considered neighbourhood, it means 
that level-set function describing the front starts to be smeared by the diffusion. Escaped particles are 
then employed to recover the mass loss, i.e., to locally reinitialize the level-set function. 

The Smoothed-Particle-Hydrodynamics (SPH) method 10.15) (see also the Godunov-SPH 12p9 ), is an 
efficient computational method proposed to give some regularity to a purely Lagrangian description of a 
phenomenon. This is done introducing suitable convolution kernels with compact supports and thinking 
each Lagrangian marker as a time-dependent Dirac measure. It follows that physical macroscopic 
quantities of interest can be obtained by simply summing the contribution of the sole particles within 
the cut-off radius of the kernel. Moreover, the spatial derivative of a quantity is easily computed by 
linearity using the gradient of the kernel. 


Finally, let us also mention the classical (W)ENO schemes [2[ lI]|l4j|2C)| (see also 7, Ch. 3]), which provide 
high-order accuracy of the numerical solution and can be profitably used as main ingredients of the proposed 
scheme. 


Paper organization 

In Section [2] we introduce the blended scheme for equations ([!]) and ^, while in Section [3] we study its 
theoretical properties. In Section [4] we compute the modified equation and we focus on some particular cases. 
In Section [5] we present several numerical tests in order to enlighten the features of our approach. Finally, we 
sketch some conclusions and perspectives. 


2. Blended numerical scheme 

To avoid cumbersome notations, let us assume d = 1, results being easily generalized to any dimension. In 
addition, let us consider the coupling of only two numerical schemes. Let f l T :=!lx [0,T] be the domain in 
which the solution will be computed, with SlcEa bounded domain and T > 0 the final time. Let us introduce 
a structured grid Q in fl T and denote by {{xi, t n )}, i = 0, ..., Nc — 1 , n = 0, ..., Nt — 1 the grid nodes, by 
Nq, Nt, respectively, the number of space and time nodes, and by Aa;, At, respectively, the space and time 
discretization steps. The space grid cell is C, = \_Xi — , aq + -§■')■ 

2.1. General idea 

Let us define the projection on the grid of the exact solution u of ([l]) 

Ui ■■= ^ ^ u(x,t n )dx, i = 0, ...,N C - 1, n = 0,..., N T - 1, (4) 

and in particular the discrete initial condition 

Ui = f u( x ) dx , i = 0,...,N c -l. ( 5 ) 
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Let us also denote respectively by {V”}i >n two discrete solutions to ([!]) at ( Xi,t n ), obtained by two 

(given) explicit-in-time numerical schemes of the form 


W" +1 =Si[W n ] 

w° = u° 


v n+1 = s 2 [v n ] 
v° = u° 


n = 0,... ,N t - 1, 


( 6 ) 


where S lt S 2 : -> R Nc , W n := (W $,..., and V n := (V£,..., V^-J. The treatment of boundary 

conditions is included in the definitions of the Si s. 

We introduce now the blended scheme defining the two new discrete solutions and {U/ 1 };,™ to ([I]), 

with W n := (W 0 n ,..., W^ c _ t ) and V n := (U 0 ”,..., V#^), by coupling and S 2 as follows 


W n+1 = A<Si \W n ] + (1 - A)S 2 [U n ], 
V n+i = (i _ M )5i[W ra ] + p$ 2 [V' n ], 


n = 0,..., N t — 1 


(BS) 


for some A,/i £ [0,1] and initial conditions W° = V 0 = 17° as defined in (|5|. The coupling parameters X, p 
manage the convex combination between the schemes. Some remarks are in order: 

(i) if (A, yu.) (1,1) the coupling occurs at each time step , hence we do not get just a trivial interpolation 
between the two independent numerical approximations. On the contrary, it is quite hard to say which 
behavior will be actually observed in the end. Moreover, there arises a natural question about which 
couple (\*,p*) would lead to the best result. We will try to give a partial answer to this question; 

(ii) if (A,/z) = (1,1), the two schemes evolve independently and W n = W", V n = V”; 

(iii) if either A = 1 or p = 1, the blend is one-way and only one solution is affected by the other; 

(iv) if (A ,p) = (0,0) the new scheme alternates the two original schemes at every time step, starting from 
S 2 in the case of W n and 5| in the case of V n ; 


(v) if A € (0,1) and p = 0, (BSl reduces to a single two-step scheme 


W 1 = A<Si[[/°] + (1 — A)S 2 [f7°] 

W n+1 = ASi[W n ] + (1 - A)S 2 [Si[ lb"" 1 ]], n > 1; 


(vi) if A = 1 — p we have W n = V n and (BS I reduces to a single scheme 

W n+1 = ASi[W n ] + (l-A)S 2 [W"], n > 0; 


(vii) generalization to multi-blend is straightforward. If we have C schemes S±,... ,Sc, the blended scheme 
is 

W n+ i’(D \ / A} ... A l\ ( 


W n +UC) J y Af ... A £ J \ S c [W n ’( £ )] 

where the coupling parameters {A)} are such that A* = 1, for all* = 1,..., £; 

(viii) a more sophisticated blend could allow (A, p) to depend on space, time, and the solution itself. In this 
way we can construct ad hoc recipes wherever the solution shows some critical behaviour. 



2.2. Multiscale coupling 

In this section we transform the blended scheme ( |BS| ) in a multiscale scheme, by choosing an Eulerian 
approximation for (e.g., Upwind, Lax-Friedriclis, (W)ENO, Godunov, etc.) and a Lagrangian approximation 
for S 2 . First, we consider a number of Lagrangian particles which move in Q driven by the velocity field A. We 
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assume that, at the initial time t = 0, we have Np particles. We denote the (positions of the) particles at time 
t by V a {t), for a = 0,..., N P — 1. We also denote by N P c = ^ £ Q the number of particles per cell. 

Initialization 

At time t = 0, the particles are uniformly placed in 11, i.e., 


V a (0) = x 0 + aAp, a = 0,...,N P -l, Ap=-^-^—^-. (7) 

Moreover, each particle is endowed with a “weight”, i.e. a portion of the mass transported by ([!]). Denoting by 
Ad a (f) the mass assigned to the cuth particle at time t, we set 

M a (0)=u(V a (0))^. (8) 

Note that, in a purely Lagrangian framework, there is no need to consider a time-dependent mass for each 
particle because it remains constant during the evolution. Therefore, one can initially distribute the particles 
in space according to u, interpreting the initial condition as the probability density function of the mass. Here 
instead we define the mass Ada as a time-dependent variable because of the combination with the Eulerian part 
of the algorithm, which will modify the masses at each time step. 

Note that one can avoid to put particles in regions in which u = 0, thus saving CPU time. For simplicity, we 
put particles all along fl, regardless of the values of u. 

Particles evolution 

As already mentioned in the introduction, the dynamics of each particle a is given by the ODE 


V a {t) = A(P a {t)). 


(9) 


( 10 ) 


It is then natural to solve ([9]) by means of an ODE solver. Focusing on one-step ODE solvers, the general 
algorithm has the form 

( VZ+ 1 = SodeCP”), n = 0,..., Nt — 1, 

l K=V a ( 0), 

where V ,£ :=V a (t n ). 

Computation of the Lagrangian density 

In order to combine the approximate Eulerian density with the Lagrangian computation we need a suitable 
averaging of the masses of the particles TVs, so to recover an approximate Lagrangian density. To this end we 
set 

V^:={S 2 [P n+ \M n }) i = ^ x Y, M «- (H) 

a : VZ +1 eCi 

where Ad™ = Ai a (t n )- Let us explain this definition: at any time t n and for any cell Ci, we find the particles 
falling in the cell Ci, summing all their masses and dividing by the area of the cell. Note that V n+1 is computed 
by the updated positions of the particles V n+1 and the last available values for the masses. 

Blending and mass update 

Given W™ +1 := (Si\W n ])i and V )” +1 as in (|TT|) at any space and time node, we are ready to run the scheme 


(BS) and compute the blended approximations W and V. However, modifying the Lagrangian density (from 
V to V) leaves the particles unchanged (both masses Ad and positions V). Since the Lagrangian density is 
computed at each time step by means of Ad and V, see we are going to lose the coupling effect after every 
time step. To fix this, we adopt the following strategy: at every time step, we modify the masses Ad™ according 
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to the modified Lagrangian density V n+l . More, precisely, we denote by A" the number of particles falling in 
the cell Ci at time t n and update their mass by setting: 


M n+1 

J 1 OL 


= M n „ + 


^(v i n+ 1 -(S 2 [V n+ \M n ]) l ), ifA?>0 
0, otherwise. 


( 12 ) 


Obviously, since the Lagrangian density is constant on each cell, all the particles falling in the same cell are 
equally affected. Note that, as opposed to the masses, the positions V a are not affected by the coupling. 


Final algorithm 

Summarizing, the complete algorithm for one time step (W n ,V n ; V n ,Ai n ) 
is the following: 

(1) Compute Wf +1 = (5i[W"])i. 

(2) Compute V™ +1 by @ (all a’s). 

(3) Compute F” +1 = (S 2 [V n+1 , M n ]) t by ||n) (all i’s). 

(4) Compute W n+1 and V n+1 by ( |BS| ). 

(5) Compute A4” +1 by (12) (all a’s). 


( w n+1 , V n+1 ; V n+1 , A4” +1 ) 


Remark 2.1. (Partial coupling) If p, = 1 the algorithm simplifies since there is no need to correct the masses, 
see (Tt2|). 


Remark 2.2. (Parallelization) The Lagrangian part of the multiscale algorithm can be easily and efficiently 
parallelized on both distributed and shared memory architectures since particles do not need to communicate 
with each other. In particular it is not needed to find the neighboring particles of each particle (as required in, 
e.g., SPH method). 

2.3. Extension to one-dimensional nonlinear conservation laws 

In this section we show how the blended scheme can be extended to equation ([3]). The Eulerian-Eulerian 
coupling does not need any modification and then it is no further discussed. The multiscale coupling (Section 
2.2) needs instead some corrections. In conservation laws, the velocity A of the particles (not that of the 
characteristic curves) can be recovered from the flux / as A(u) := (if a / 0) and equation © must be 
replaced by 

V a {t) = A(u(V a {t),t)). 

This means that the Lagrangian scheme needs the density u to update the positions of the particles. The 
blended scheme provides at each time step two approximate densities, namely W n and V n . Hence we solve 
either 

Pa(t n ) = A(W£"J or p a (t n ) = A(V£'J , (13) 

where j a ,n is the cell where the a-th particle falls in at time t n . Note that, choosing W n , we introduces 
an additional coupling: the Lagrangian scheme now depends on the Eulerian scheme even for an uncoupled 
evolution (A ,/j.) = (1,1). 


2.4. Estimation of coupling parameters (A, fi) by Richardson extrapolation 

The blended scheme ( |BS[ ) is actually a family of numerical schemes indexed by two coupling parameters 
(A,/x) £ [0, l] 2 - Generally speaking, we expect that only some couples lead to a solution W or V which is 
better than both solutions W and V obtained by the uncoupled schemes. The question arises how to find an 
advantageous couple a priori , before running the computation and without the knowledge of the exact solution. 
In some simple cases one can compute the optimal parameters analytically, but in general this is impossible. In 
order to provide a general method, we adopt the following strategy. 
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First of all, let us denote the solutions of the blended scheme (BSl by {W"(A, [i)}i,n and {F"(A, n)}i, n in 
order to make it explicit the dependence on the coupling parameters. Moreover, let us define the L 1 errors at 
final time as 


E l [W{ A,m)] “ E \ W ? T ^V) - U? T Ax, E 1 [V(X,fi)] := E \,») - U t 


N'T 


Ax, 


(14) 


and define the two optimal L 1 couplings as 


(A *,H*)w := argmin E 1 [W (A, /z)], (A*,/z*)y ■= argmin£' 1 [F(A,/x)]. 


(15) 


Following the classical idea of Richardson extrapolation, we can run the blended scheme on two (coarse) grids, 
Q’ and Q ", thus obtaining, for each (A, n), four solutions W(X, ii;Q'), V(A W(X, and V(X, 

Then we define the two Richardson L 1 error indicators as 


6*(X,Lr,G',g") := E | W? T (A, /x; Q') - W^(A,/z; Q") 

i&Q' 

<^u(A>M; S',Q") := E I^CA,^;^)-V^5(A ^G") 

ieS' 


Ax, 

Ax, 


(16) 


where h(i) gives the best correspondence from node i in Q' to a node in Q". The estimates for (A *,n*) are then 
given by 


( X R ,n R ) w ,g',g" ■= arg min S^(X, /z; Q', Q"), 

A ,fi 

(A R ,H R )v,g',g" := argmin5y(X,n;G',G"). 


(17) 


Remark 2.3. It is possible to avoid an exhaustive search in [0, l] 2 in (JTt]) by employing a descend method in 
the space (A, /z) starting from some initial guess (Ao, fio)- One can also progressively refine the grid (in the space 
(A , fj,)) around the argmin found in the previous level. 

Remark 2.4. A necessary condition for the method to work is that (A*, fi*) is stable with respect to the grid 
size. In Section [5] we will investigate this property. Moreover, we expect the method to be sufficiently fast if 
performed on very coarse grids, in particular much coarser than the final grid where the actual computation 
must be performed. 


3. Theoretical analysis 


We address in this section the main theoretical issues for the blended scheme solving equation ([!]) with d = 1. 
Note that, rather than perform a complete convergence analysis (which apparently requires more technical 
subtleties), we will analyze the scheme in the two situations which seem to be more effective and that are 
actually used in numerical examples, i.e., the coupling of two Eulerian schemes and the use of a Lagrangian 
scheme to correct an Eulerian scheme in a multi-scale approach (i.e. A € [0,1) and /i = 1). 


3.1. Mass conservation 


We first prove that the scheme (BS) is conservative, provided the two constitutive schemes are so. Let us 
use again the notation 


= {S 1 [W n ]) i , v? +1 = (S 2 [V n ])i (18) 

for a single time step of the uncoupled schemes (with a slight abuse of notation in the case of multiscale 
algorithm). 
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Proposition 3.1. Assume that the solvers © are conservative, i.e. 

E wr +1 = E w - r. E ^” +1 - E ^ 


and that 


Ewf = ek° = E^- 


(19) 

( 20 ) 


Then, the scheme (BSl is also conservative. 


Proof. By definition of the scheme (BS) and (18), we have 


' E< wr 1 = A Ez w ? +1 + (1 - A) Ez v ? +1 

Ei ^” +1 = (1 - T) Ei w” +1 + t Ei ^ n+1 , 


( 21 ) 


so that, using the discrete conservation of both schemes (19), we have 


'Ei w? +1 = A Ez Wf + (1 - A) Ez 

,Ei ^ +1 = (i-M)Ei^r + MEi^ 


which implies that 


Er.EEf 


E^’E^” ’ max E^-E^ 


The proof is concluded by induction noting that (20) holds. 


□ 


Remark 3.2. The previous result remains true if both parameters A and /x depend on the step index n, but 
not if they depend on the cell index i. 

3.2. Coupling between two Eulerian schemes 

Hereafter we will assume for simplicity that the schemes iSi,<S 2 are linear, so that we can write 

W n+1 = S'iW” , V” +1 = S 2 V n , 

where ,S’-|, S- 2 are two Nq x Nq matrices. Let us also assume that each of the schemes is stable, and in particular 
that they satisfy 

Il'S'illpi < 1 + CiAt , \\S 2 \\ P2 < 1 + C 2 At, (22) 

with constants C\, C 2 independent of Ax and At, and for some (possibly different) suitable norms || • || p . (i = 1, 2) 
which could be required for the convergence analysis of the two schemes (e.g., oo-norm for a monotone scheme 
and 2-norm for a high order scheme). We also assume to use scaled l p norms defined by 


i/p 


\\U\\ P = Ax 1 / p E 


(23) 


Now, the scheme (BSl can be rewritten in matrix form as 


W n+1 \ (XI (1 - A) A (S x 0 
V n + X Ul -ti)I ,il 0 S 2 


W* 

yn 


= QS 


W ’ 

yn 


(24) 
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where I denotes the Nc x Nc identity matrix, and 

, XI (1 - A)/ 
1 (1 - H)I Hi 


s = 


Si 0 
0 S 2 


In order to analyse convergence for the scheme (BS), we define a suitable norm for the numerical solution 
(W n , V n ) in the product space M^ 17 x M. Nc as 


\\(W,V)\\ B := max(||W|| Pl , 


,)• 


Note that 


|0||b = 1 , H-S'Hb < 1 + max{Ci,C , 2 }At, 

b, since 


(25) 


which implies that (BS) is stable in the norm 

lies’ll < ||0 ||b||<S , ||b < l + max^C^Af. 


It is also immediate to verify that (BS) is consistent, and therefore convergent. In fact, assume to work in the 
norm || • ||b and apply the usual definition of consistency for a discrete solution at time t n defined by 




where U n is defined in Q for a smooth u(x,t). Taken into account that the discretization parameters are At, 
Ax, the scheme is consistent if and only if 


(U n+1 \ 

\U n+1 ) 



T At 


f &i{Ax, At)\ 
\a 2 (Ax,At)J 


(26) 


with <ti,( 72 -A 0 for At, Ax —> 0. Now, if the two schemes have respectively truncation errors t\ and r 2 , it is 
immediate to compute and cr 2 as 


(ai{Ax,At)\ _ (t!(A x, At)\ 

\a 2 (Ax,At)J \T 2 {Ax,At)j 


which results in consistency for the scheme (BS) if the two elementary blocks are consistent. We can therefore 
conclude with the following 


Theorem 3.3. Let Si, S 2 satisfy (22), and , for i = 1,2, 


Ti{Ax, At) —> 0 (Ax,At—>0). 


Then, the blended scheme (BS) is convergent in the norm || • ||b for any A, yi £ [0,1]. 


3.3. Convergence for the Lagrangian correction 

We recall that in this case we are assuming that the Lagrangian part of the scheme provides a correction for 
the Eulerian part. This requires to set /i = 1, whereas A G [0,1], and in this case we clearly have V n = V”. In 
order to have a more explicit notation, we will denote here the two components of the error at the n-th time 
step for the Eulerian and the Lagrangian part as ef, and ef, respectively. 

We start by computing the error for the pure Lagrangian part of the scheme, then we turn to the coupled 
scheme. 
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Figure 1. The characteristic map $. 


Step 1. Bound on := U n — V n 

We denote by $ : (xcrO) i —> (x,t n ) the mapping between the foot xq of the exact characteristic curve at time 
t = 0 reaching the point x at time t n (see Fig. [lj. 

Similarly, we denote by the mapping corresponding to the discrete characteristic curve, given by the ODE 


solver (10). We will make the standing assumption that this is a monotonic and invertible map (i.e., that no 


crossing between discrete characteristics occurs). In usual ODE solvers, this happens for At sufficiently small. 
We also assume that the error of ODE scheme is of order 0(At 7 ), for some 7 > 1. 

Then mass conservation and it(x, 0) = u(x) imply 


UP = 


1 

Ax 


f u{x,t n )dx = [ u{y)dy=[ u($ ^x))^ 1 (x))'dx, 

JCi A x i(Ci) Ax JCi 


where the last equality follows applying the change of variable y = $ 1 (x). 
We also define 

UP = J c u(<S>A 1 (x))(<I>pp(x)ydx 

to split the estimate as \Up — V™\ < \Up — V™\ + \Up — Up\. 

Let us recall the definition of Vp as in 


V n = V r 


1 

Ax 


E M 


where the masses are computed by backtracking the Lagrangian particles up to the initial time. Let N be the 
number of particles Pp falling in C”i, indexed for simplicity by k = 1 Note that N = O(Npc) by the 

Gronwall’s lemma. Since Pp is the evolution at time t n of the corresponding at time zero, we have 


V n 

l 


1 

Ax 


N 


E«(^) 


fc=1 


Ax 

Npc 
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Using the discrete map we can decompose the cell Ci as the union of the pre-images of uniform intervals 
centered at the particles, namely 


N 


C = U f I>A 


k =1 


Pk- 


Ax 
2 Npc 


’ Pfe ° + wt) )u ( [p°- 1. p fe ° + i)) ■ 


Applying the change of variable y = $ A 1 (x) we obtain 


N 


iU-ui<^E 


k =1 


r ” fc+ i At 

/ + 2 u(y)dy-u(P°)-^ 

Ipo JSpc 

k-i 


1 ^ Ax 


Aa; Npc 


u ( P 0 k +O 


( Ax \ 


~u(P°) 


V Npc J, 

where the last step follows from the mean value theorem and the fact that P° + i — P^_ i = Then we get 

(27) 


i uk-vn<^pL t o(^-)<c^-, 


where L. a is the Lipschitz constant of the initial datum. 
Now, by the change of variable $ A (x) = $ -1 (y) we get 


\u?-u?\< 


1 

Ax 


■u($ (y))($ ( y))'dy - / u($ (*))($ (x)) dx 

JCi 


where Ci = <I>(<I> A ( Ci )). In practice (see Fig.|2j), we trace back each point x £ Ci along the discrete characteristic 
up to < & A 1 (x), then we move forward along the exact characteristic up to y = $($ A 1 (x)) € Ci. 


t 



Figure 2. Combining exact and discrete characteristic maps $ and <F A . 

Since x and y share the same foot, the error estimate for the ODE solver provides a function (j> such that 

y = x + <j>{x)0{AC ), 

namely a mapping between C t and Ci, which in turn gives 


W-u? \< 


Ax 


Ci 


u(<E> 1 (x + cj)(x)0(At 1 )))(<& 1 (x + <t>(x)0{AC))Y — u($ 1 (x))(<l> 1 (x)) / dx + O(AC) < 
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< C^~ [ \d>(x)\0(AF) + 0( AC) < CAP . 

A x JCi 

Combining (27) with (28), for any norm in the form (|23|), the following bound holds: 


\el\\ p <C(At'y + 


Ax \ 


Npc 


) 


(28) 


Step 2. Bound on el := U n — W n 

Writing (BS) between n — 1 and n. and substituting the second relationship in the first one, we get 

W n = XSiW 71 - 1 + (1 - A) ( U n - el ), 
which gives in turn, for the error on the component W n , 


C E 


= U n — W n = \[U n — S'llf” -1 ] + (1 — A)e2 = 

= A {(U n - S it/”” 1 ) + (Sit/ 71 " 1 - S'iW 71 " 1 ) ] + (1 - X)e n L . 


(29) 


Using now the consistency error t\ of the Eulerian part of the scheme, we can write the first term in square 
brackets as 

U n _ SiU n -1 = Atn(Ax, At) 

which gives, when used in 


e n E = X [At n(Ax, At) + (S^ET*" 1 - W n ~ x )) ] + (1 - X)e n L . 
Passing to the norms, we obtain therefore 

+ AAf ||7i(Ax, At)|| p + (1 — A) ||e2|| p • 


\ e E lip < A(1 + CiAt) ||e 


n— 11 
E I 


Defining now ( := A(1 + C\ At) and iterating back the estimate, we obtain 


\ e E lip < C|| e E 


n—1 1 


< c 2 h r ! 


i + AAt||r 1 (Ax,At)|| p + (l-A)||c2|| p < 

+ AAf ||ri(Ax, At)|| (1 + C) + (1 - A) (||e£|| + C 


"p 
n—1 


< • • • < 


< C" II4IL + AAf IlntAx, At)|| p J2( k + (l-X)J2C 


k || n—k II <r - 

i-i ||p — 


k—0 


k—0 


< C \Ue L + AAf ||Ti(As, At) || + (1-A)max 


n—k | 


n—1 


£c‘- 


k =0 


Taking into account the definition of £, the last term in (33) can be bounded as 

OO 

if A < 1 


n—1 




k =0 


1-(A + 6) 


k =0 


e c ' T - 1 
Ci At 


(30) 

(31) 


(32) 


(33) 


if A = 1, 

where the first case clearly holds for At small enough to have A(1 + CiAt) < A + <5 < 1. Dropping the trivial 
case A = 1 (which gives two uncoupled schemes in (BS)), we obtain therefore, for A < 1 and At small enough, 
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the estimate 


4IL < ||4|L + C'(A) ( At ||ti(A.t, At)||„ + max| 


n—k I 




Step 3. Conclusions 

Finally, assuming that the initial condition has exact cell averages (i.e., that 11 11 D = 0), and recalling that 

Npc = = !VpO(A:r), we obtain that the scheme converges and the error e n = ( 


Nc 

| b, the error bound 


1 p 

e u> £ l) satisfies, in the norm 


<C(At\\T 1 (Ax,At)\\+AC + 


Ax 


•1 = 


C( A(||n(Ai,A()|| p + AC + — ), (34) 


Npc J 

for any n > 0 such that t n £ [0, X]. We can then summarize this analysis in the following 


Theorem 3.4. Let Si satisfy (22), and assume that ti(Ax, At) —»• 0 for Ax, At —»• 0. 

Then, for any A £ [0,1], the blended scheme (BS), as detailed in Subsection 2.2 is convergent in the norm 
|| • ||b as Ax, At -a 0, Np —► +oo, and the estimate (34) holds. 


4. Modified equations 

In this section we derive a system of modified equations of |l]) (see 13" Sect. 11.1]) associated to the scheme 


(BS), in order to give an interpretation of the effect of the coupling between the schemes from a continuous 
point of view. 


4.1. The general case 

For our purposes, it is convenient to write first the scheme in the following implicit form, obtained by 
straightforward algebraic manipulations in the case X ^ 1 — p: 


W n+1 -5 1 |W”1 = —— —<V n+1 - W n+1 ) 
X + p-1 


V n+ ^S 2 [V n ] = 


1 - n 
A + fi — 1 


u = 0 ,..., Nt — 1 • 


(35) 


(IF n+i - V n+i ) 


Now, let us assume A(x) = A and that, for a given Ax/At, the schemes S\ and S 2 have respectively the modified 
equations 

d qi w d q2 v 

w t + Aw x = viAx Pl — -, v t + Av x = v 2 Ax v2 


dx 9 1 ’ 


drS 2 ’ 


for some p\,q\,p 2 ,q 2 . Dividing by At in (351, using a Taylor expansion of the right hand sides up to first order 
in At and restricting to the leading terms, we end up with the following system for w and v: 


w t + Aw x = V\ Ax p 
v t + Av x = v 2 Ax P2 


i d qi w 
dx q i 
d q2 v 


1 — A 


At A -(- p — 1 

1 1 -n 


(v — w) + 


1-A 
A + /i — 1 


(vt - w t ) 


dx q2 At X + pL—1 




(36) 


w{x, 0) = v(x, 0) = u(x) 

We clearly see that the dominant terms involving 1/At act as reaction terms, depending on the sign of the 
differences ±(v — w). We expect that, at the very beginning, w and v start to deviate from each other. This 
immediately activates both the reactions, driven by the weights 
solutions close to each other. 


1 —A 
— 1 


an d respectively, which push the 
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Remark 4.1. The reaction terms appearing in the modified equation 36 are similar to those appearing in 16 


Eq.(10)], where the authors combine the Lagrangian SPH solution with a Finite Volume scheme. There, a 
“chimera”-like approach is used to force one solution towards the other, in order to minimize discontinuities 
and artifacts. 


4.2. The particular case A = 1 — p 

If A = 1 — /x we have W n = V n (cf. (vi) in Section [2]), and with some minor manipulation we can rewrite 
(BS) as 

W n+ 1 - w n _ Si [W n ] - W n S 2 [W n ] - w n 

— A-—-b (1 — A) - 


At At v ' At 

Restricting again to the leading terms, this yields the modified equation 


f) 

w t + Aw x = ( Aui Ax Pl a + (1 — A)u 2 Ax p " 


dx 91 


<9x® 


which shows a final dispersive term obtained by a mere convex combination of the two original ones. 

Note that, if the two schemes are not of the same order, i.e., if p\ ^ p 2 , then the scheme of lower order 
becomes asymptotically dominant if A and p are kept constant along a refinement. On the other hand, if a 
nontrivial optimal combination exists, then both A and p are expected to depend on Ax. 


Example 1: Lax-Wendroff + Beam-Warming 

Assume that A = 1 — p, and that we are coupling the Lax-Wendroff (LW) and the Beam-Warming (BW) 
schemes 13 Chapt. 10] in the simple case A(x) = A > 0. As a consequence of having dispersive terms of 


the same order, but different sign, these two schemes are known to have opposite behaviour with respect to 
numerical dispersion. In fact, the first one causes oscillations behind a discontinuity, while the second one ahead 
of it. We have therefore some chance to compensate the two opposite responses. 

Denoting the Courant number by /3 = A At/ Ax, we have for the LW scheme the modified equation 


while the BW scheme gives 


AAx 2 , 

w t + Aw x =--— (1-p ] w xxx 


AAx 2 

w t + Aw x =--— (2 - 3/3 + /3 2 ) w x 


Now, the coupled scheme has a modified equation given by 


A Ar 2 

w t + Aw x = —[A (l - /3 2 ) - (1 - A) (2 - 3/3 + /3 2 )] w xxx , 


and hence, imposing the term in square brackets to vanish, we obtain the condition 3A(1 — /?) = /3 2 — 3/3 + 2, 
which finally gives the value of the optimal parameter A* = Clearly, this procedure makes the leading 

dispersive term vanish, thus the modified equation of the blended scheme will now have a nonzero term in 
w xxxx . A straightforward computation shows that the resulting scheme is nothing but a third-order Upwind (or 
semi-Lagrangian) scheme obtained by interpolating the numerical solution W n at the foot of a characteristic 
fine via a symmetric Lagrange polynomial of third degree. 


4.3. The particular case V = U 

Some additional considerations can be done in a even more special case, namely assuming that the scheme 
<5>2 corresponds to the exact solution (i.e. V = U): 


V," = u(xi, nAt) = u(xi — An At). 
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This analysis is interesting because gives some insights about the coupling with nondiffusive schemes. In 
particular, it shows the limit behavior of the coupling with a Lagrangian scheme when the number of particles 
tends to infinity. 

Note that the blending of V with 6>i does not make sense in this case, being V exact, therefore we keep /i fixed 
to 1. Again, due to the numerical dispersion of Si, one could expect that the blended solution W still shows 
a diffusion, just mitigated by the exact solution. Surprisingly, this is not what is observed in the numerical 
simulation, which instead shows that the numerical diffusion disappears after a short transient. This behavior 
can be explained by looking again at the modified equations. After a rearrangement of the terms containing 
Wt, the system (36) reduces to 


f »( + A Aw x = XvAx p ^^ + u - w) + (1 - A )u t 

^ ut + Au x = 0 (37) 

l. w(x , 0 ) = u(x , 0 ) = u(x) 


for v = Vi, p = pi and q = qi- It is easy to verify that the difference r := w — u satisfies 


1 - A d q d q 

r t {x,t) + A Ar x (x,t) = --^-r(M) + XuAxP Q xq r ( x A) + XvAx p —u(x,t), 


(38) 


with the initial condition r(a;,0) = 0. Performing a Fourier transform with respect to the space variable, we 
get: 

r(£, t) = -A Ai£r(£, t) - t) + AuA x p (i£) q r(£,, t ) + XvAx p (i£) q u(£, t ) = 

= (~XAi£- ^^ + XuAx p (i^f(^t) + XuAx p (i^ q u(^t), (39) 


where £ is the Fourier variable and f denotes the Fourier transform of r. We observe now that (39) represents 
a family of decoupled linear ODEs indexed by £ with source terms given by XvAx p {i^) q u{£ t ,i), and that the 
coefficient in brackets has a negative real part for any £, provided the dispersion operator is dissipative (as it is 
usually the case in stable schemes). Then, we infer that, for t —> oo, the function r(£, t) (and hence, the function 
r(x,t )) remains bounded in L 2 provided has itself a bounded 2-norm. This shows that the numerical 

diffusion of the blended scheme is eventually balanced by the reaction term, and no further deteriorates the 
solution. This observation is confirmed by the following numerical example. 


Example 2: Upwind + Exact 

Let us couple the Upwind scheme (Si) with an exact scheme (S 2 ) with U = [0, 20], T = 10, A = 1, Nq = 300, 
Nt = 800, and initial condition 


s f \ (1 + cos(tt(:e - 2))) , if I* — 2| < 1, 
u(x) = < ^ (40) 

[ 0, otherwise. 

Here the Courant number is equal to 0.1875. Fig. [3](a) shows the approximate solution W at final time for 
different values of A. For A = 1 (pure UPW) it is visible the well-known diffusive behaviour of the scheme. 
Switching on the coupling, namely assuming A < 1, the behaviour changes completely. After a short transient 
during which the solution is smeared out, the solution is frozen and propagates at correct speed forever, with 
no further change in the shape. Fig. [3j]b,c) shows the evolution of the maximum of the solution and the L 1 
error for different values of A. For A = 1, the maximum decreases, as a consequence of the numerical diffusion. 
For A < 1, the maximum decreases for a short time, and then it comes to a regime state. Analogously, the L 1 
error is bounded. 
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(a) (b) (c) 

Figure 3. Toy test II: (a) the exact solution u and the approximate solution W at final time, 
for /i = 1 and A = 0.997, 0.99,1. (b) max and (c) L 1 error of the approximate solution IT as a 
function of time, for /i = 1 and A = 0.8, 0.95, 0.99, 0.995,1. 


5. Numerical tests 


In this section we present some numerical tests, showing the behavior of the blended scheme (BSl and 
confirming the theoretical results described in the previous sections. We always assume that the CFL condition 
applies, i.e. the time step is chosen in such a way that 


At = 


, Ax 

'PIu 


(eq. 


At = /3 


Ax 

ii/JU 


(eq. 0), 


(41) 


where /3 < 1 is chosen small enough to guarantee the stability of the scheme in use. Recalling definition (14), 
in the case of Eulerian-Eulerian coupling we define the reference error E^ e{ as 


E] ei := min{.E 1 [W(l, 1)], S 1 [F(1,1)]} , 


(42) 


i.e. the minimum error one can achieve with the original uncoupled schemes. We also define the following subset 
of the parameter space 


*[ W ] ■■= {(A,M) € [0,1] 2 : E'[W{\,n)} < E^}, 
$[F] := {(A,/z) € [0, l] 2 : E l [V{X, (i)] < E^}, 


i.e. the couples (A, p.) which improve the approximation with respect to both the uncoupled schemes. 

In the case of multiscale coupling we always set /x = 1 (partial coupling) since it seems more convenient, both 
from the accuracy and CPU time point of view. Moreover, to be fair we define 

E] ei := E\W(l, 1)] and $[W] := {A e [0,1] : E 1 [W(X, 1)] < E^ { }, (44) 


using only the Eulerian density W(l, 1) and not V(l, 1). This choice is motivated by the fact that the recon¬ 
struction of the Lagrangian density proposed in (11), despite it is a natural choice in the framework of the 
blended scheme, it is not in general the best approximation one can achieve (see, e.g., the SPH method). As a 
consequence, it is quite easy to get more accurate approximation of the Lagrangian density. 

In order to use Richardson extrapolation, we introduce a scale parameter s € (0, |] to coarsen the grid in 
a proper way. Given parameters Nc, Nt for a reference grid Q, we define (denoting by [■] the upper integer 
part) the grid Q' such that N' c = [s_/Vc], Nj, = [sJVr], and the grid Q" such that Nq = 2[siVc] = 2 N' c , 
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Nfp = 2 [sNt] = 2 Np. We also define Np = [slVp] and Np = 2[slVp] = 2Np to scale the number of particles, see 
section [3] 

Hereafter we use the following acronyms: UPW = Upwind (I order), RLW = Richtmyer two-step Lax- 
Wendroff (II order) 13j Sect. 12.2], WEN02 = Weighted Essentially Non-Oscillatory with linear interpolation 
and TVD Runge-Kutta 3 approximation in time (small constant e = 1CU 9 ) (III order) 

Forward Euler (I order). 


14 


EE = Explicit 


5.1. Test 1 


Test 

Si 

s 2 

u(x) 

n 

T 

A(x) 

N c 

N t 

P 

1 

RLW 

UPW 


[0,20] 

2.3 

X 

1200 

3000 

0.92 


In this test we blend RLW and UPW in case of the positive space-dependent velocity A(x) = x. The exact 
solution is given by u(x,t ) = u{xe~ t )e~ t . Fig. [ 4 ] shows (a) the level sets of the function (A,/r) —> E 1 [W(A, /i)] 
and (b) the region $[W]. In this case one can note that the level sets of the error are straight lines and one can 




Figure 4. Test 1: (a) level sets of the function (A ,fi) —► E 1 [W(X, n)]. (b) Region <&[W] 
(white). The red square indicates the global minimum point (A*,/U), the green circles indicates 
the global minimum points (A*,//*) on coarser grids (s = 2 ’ 3 > 4 > 5 )’ an d magenta cross 
indicates the Richardson minimum point (X R ,/i R ). 


obtain the same solution varying the coupling parameters while keeping their ratio fixed. This implies that the 
line A = 1 — /z nearly cuts all the levels sets, and, therefore, restriction to such a particular case (much simpler 
to study analytically since the system is reduced to a single equation) does not lead to a loss of generality. 

The reference error E^ cf in this case is equal to 0.1463, while the minimum error that can be achieved is 0.0816 
(-44.22%), corresponding to the couple (A*, /U) = (0.8533, 0) computed by an exhaustive search in the parameter 
space. Richardson extrapolation method (with s = |) suggests instead the couple (X R ,n R ) = (0.29,0.1), which 
is fully inside <&[W] and leads to the error 0.117 (-20.03%). Fig. |4](b) also shows the best coupling parameters 
(A*,/U) computed coarsening the grid with s = 5 , 4 , g- It can be seen that the best coupling parameters do 

not depend very much on the grid and move monotonically with respect to the scaling from right to left. 

Fig. [5]shows the exact solution U Nt together with (a) the two uncoupled solutions W Nt ( 1,1) and V Nt (1, 1), 
(b) the best coupling W Nt ( A*,/U), and (c) the best Richardson coupling W Nt {X r , fx R ) at final time. Note that, 
although W Nt (X r , ^l r ) is worst than W Nt ( A*,//*) in terms of L 1 error, it does not show spurious oscillations. 
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Figure 5. Test 1: comparison with the exact solution, (a) 1) and V Nt (1, 1). (b) 

W Nt {\*, n*). (c) W Nt (\ r ,h r ). 


This is due to the fact that oscillations are sensitive to the grid and then Richardson method easily discards 
such a solutions. 

To better quantify the improvement of the blended scheme with respect to the uncoupled schemes, we found 
the grid needed by the i 1 -best uncoupled scheme (RLW) to achieve the same accuracy of W Nt (X r , fi R ) and 
W Nt (A*, /li*). In the first case we need Nc = 1704, Nt = 4260, while in the second case Nc = 3240, Nt = 8100, 
corresponding respectively to a refinement of the grid of a factor 1.42 and 2.71. 

5.2. Test 2 


Test 

5, 

s 2 

u(x) 

n 

T 

A(x) 

N c 

N t 

Np 

P 

2 

UPW 

EE 

X \bi]W 

[0,20] 

2.3 

X 

1200 

3000 

5 x Nc 

0.92 


In this test we run the multiscale version of the blended scheme coupling UPW and EE in the same setting 
of Test 1. The UPW error, which will be taken as the reference one, is equal to 0.1771. The minimum error 
that can be achieved is 0.0204 (-88.48%), corresponding to A* = 0.992 computed by an exhaustive search in the 
parameter space [0,1] x {1}. Fig. |6| shows the function A — > E 1 [W(X, 1)] cut by the reference error. Richardson 
extrapolation method (with s = |) suggests X R = 0.99 which leads to the error 0.0208 (-88.26%). 

Fig. [ 7 ] shows the exact solution U Nt together with (a) the two uncoupled solutions W Nt (1, 1) and V Nt (1, 1), 
and (b) the best Richardson coupling W Nt (X r , 1) at final time. It can be seen that EE, as opposed to UPW, 
is rather good in terms of diffusion and in capturing discontinuities, but it has several oscillations. The size 
of the oscillations depends on Np, while their number depends on Nc- The vector field A(x) = x tends to 
separate the particles (but in this case they remains equispaced), therefore the oscillations become larger as 
the time goes by. Note the the Lagrangian density V(l, 1) is 0 in those cells not covered by any particle. The 
blended solution W Nt (X r , 1) combines nicely the advantages of the two schemes, showing little diffusion and 
tiny oscillations (which have a different nature with respect to those of high-order Eulerian schemes). For the 
sake of comparison, we report that WEN02 scheme reaches the same accuracy of UPW+EE on a grid refined 
by a factor 4.2. 

We also tested the possibility of localizing the particles around critical zones, namely the discontinuities. At 
initial time, we locate a few particles across the right discontinuity, and then we activate the blend (with A = 0) 
only in the cells which are between the first and the last particle (more in general, inside the support of the 
Lagrangian density). This is needed in order to distinguish the two kind of no-particle zones, the one due to 
the particle rarefaction and the one uncovered by user’s choice. In Fig. [8](a) we show the result of the blended 
scheme in this case and the region actually covered by particles. The right discontinuity is well caught despite 
the small number of particles. 
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Figure 6. Test 2: function A —» E 1 [W(X, 1)] cut by the reference error, (a) A € [0.6,1]. (b) 
A £ [0.89,1], The red square indicates the global minimum point A*, the green circles indicates 
the global minimum points A* on coarser grids (s = \ , g, g), and the magenta cross indicates 

the Richardson minimum point \ R . 
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(b) 


Figure 7. Test 2: comparison with the exact solution, (a) W Nt ( 1,1) and V Nt (1, 1). (b) 
W Nt (X r , 1). Solutions are downsampled for a better viewing. 


For the sake of completeness, we also report the outcome of the blended scheme whenever the Lagrangian 
solution is regularized by the Eulerian one. In Fig. [8][b) we show the solution V obtained with (A, /x) = (1, 0.3). 
In this case the blend puts together, rather than cancels, the drawbacks of the two uncoupled schemes, leading 
to a solution which is both diffusive and oscillating. 
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Figure 8. Test 2: (a) W Nr { 0,1) with localized particles. Green circles are the particles, (b) V Nt (1,0.3). 
5.3. Test 3 


Test 

Si 

s 2 

u(x) 

Q 

T 

A(x) 

N c 

N t 

TVp 

13 

3 

UPW 

EE 

sin(e 2 720 ) 

[0,7T] 

1 

sin(x) 

600 

200 

1 x Nc 
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In this test we run the multiscale version of the blended scheme coupling UPW and EE in case of a positive 
space-dependent velocity A(x) = sin(x) (for x € O = [0,7r]). The exact solution is 


u(x, t) = u(j(x, t)) ^ tan (l + e 2t ) sin(y (x,t)) + e 4 cos (7 j(x, t) = 2 arctan 4 tan^^. 

The reference error (UPW) in this case is equal to 0.2591. Fig. [ 9 ] shows the exact solution U N ' r together with 
the two uncoupled solutions W Nt ( 1,1) and V N ' r (1,1) at final time. In this case the particles move rightward 
accumulating on the right side of the domain. The minimum error that can be achieved is 0.0731 (-71.79%), 
corresponding to A* = 0.93. Richardson extrapolation method (with s = |) suggests instead \ R = 0.916, which 
leads to the error 0.0742 (-71.36%). 

Fig. [To] shows four blended solutions obtained with different choices of A. For low values of the coupling 
parameter (A = 0.8) the solution is accurate where the number of particles is large (right side) but not where 
there are no or a few particles (left side). The situation is reversed for large values of the coupling parameter 
(A = 0.99). The Richardson solution A = X R = 0.916 reaches a good compromise between the two cases. On 
the other hand, an even better result can be obtained by means of a variable-in-space A = A(A) (see Section 2.2 
for the definition of A). Last plot in Fig. 10 shows the result of an hand-tuned decreasing function A(A) which 
leads to a good accuracy in both sides of the domain. 
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(b) 


Figure 9. Test 3: uncoupled schemes in [2.3,7r]. (a) W Nt {1,1). (b) V’ JVr (l, 1). 


5.4. Test 4 
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In this test we solve equation @ with f(u) = u( 1 — u) (LWR model for traffic flow). The exact solution is 
u(x,T) = — + | in [2,6] and u(x,T ) = 0 elsewhere. A rarefaction and a shock are present. The reference 

error (GODUNOV) in this case is equal to 0.0839. Richardson extrapolation method (with s = \) suggests 
X R = 0.956 which leads to the error 0.0317 (-62.22%). Here we used W to evaluate the particles’ velocity, 
shows the exact solution U Nt together with (a) the two uncoupled solutions W Nt ( 1,1) and 


see 


©■ Fig. 


11 


V' JVr (l, 1), and (b) the best Richardson coupling W Nt (X r , 1) at final time. 


Conclusions and perspectives 

We proposed a new method for coupling two (or more) explicit schemes approximating advection equations 
and scalar conservation laws. We provided some convergence analysis and we showed that effective couplings can 
be obtained in practice using a standard Richardson extrapolation. Moreover, several numerical tests confirmed 
that, whenever the constitutive schemes are respectively an Eulerian (macroscopic) and a Lagrangian (particle 
based) scheme, the blended scheme gives rather accurate results despite its Eulerian part gives only a low order 
approximation. This is a remarkable advantage, since it is well known that implementing high-order solvers is 
much easier for ODEs than for PDEs. 

It is interesting to note that the blended scheme can be put in relation with very different schemes proposed 
in the literature, see Section [l] 

• If the coupling parameters are allowed to depend on space and/or on the solution itself, we can blend 
Eulerian schemes of different orders to recover the constitutive ideas of the filtered schemes. 

• Coupling an Eulerian scheme with a pure Lagrangian scheme and localizing particles, we recover a 
scheme that resembles that of the particle level-set method. 

• Annealing the Eulerian scheme but keeping alive the coupling with (A, p) = (0,1), we recover a scheme 
that resembles that of the particle-in-cell method. 

• The multiscale blended scheme (with /i = 1) seems to follow the opposite philosophy with respect to 
the smoothed-particle hydrodynamics method. While the latter regularizes the Lagrangian solution by 
a suitable choice of smoothing kernels, the former “Dirac-izes” the Eulerian solution. Moreover, the 
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Figure 10. Test 3: W Nt ( A, 1) in [1.3, 7r] for different choices of A. from top to bottom: A = 0.8 
( E 1 = 0.1085), A = X R = 0.952 ( E 1 = 0.0764), A = 0.99 (. E 1 = 0.1517), A = A(A) (E 1 = 
0.0396). 
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Figure 11. Test 4: comparison with the exact solution, (a) W Nt (1,1) and V Nt {1, 1). (b) W Nr (X R , 1). 


computation of the optimal blending parameter resembles the computation of the optimal supports for 
the smoothing kernels. 

Finally, let us mention that one of the most promising and yet unexplored feature of the blended scheme 
is the parallelization on distributed-memory architectures. Particle evolution is indeed embarrassingly parallel 
(within a time step), since there is no need to locate the neighbouring particles of a given one (which can be a 
severe bottleneck if particles live in different memory units). 
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